function [statistical_results, fig_statistics] = statistical_channel_models(num_samples, k_factor_db)
% 统计信道模型分析
% 输入参数：
%   num_samples - 样本数量
%   k_factor_db - 莱斯因子 (dB)
% 输出参数：
%   statistical_results - 统计结果
%   fig_statistics - 图形句柄

% 添加路径
addpath('../Common');
colors = color_definitions();

%% 并行计算设置
[use_parallel, pool_info] = parallel_manager();

% 参数设置
k_factor_linear = 10^(k_factor_db/10);

%% 生成衰落样本
if use_parallel
    % 使用并行计算生成衰落样本
    fprintf('使用并行计算生成衰落样本...\n');
    
    % 预分配临时结果存储
    temp_results = cell(4, 1);
    
    parfor i = 1:4
        switch i
            case 1
                % 瑞利衰落
                rayleigh_temp = sqrt(randn(1,num_samples).^2 + randn(1,num_samples).^2);
                rayleigh_temp = rayleigh_temp / sqrt(mean(rayleigh_temp.^2)); % 归一化
                temp_results{i} = rayleigh_temp;
            case 2
                % 莱斯衰落
                rician_temp = sqrt((sqrt(k_factor_linear) + randn(1,num_samples)).^2 + randn(1,num_samples).^2);
                rician_temp = rician_temp / sqrt(mean(rician_temp.^2)); % 归一化
                temp_results{i} = rician_temp;
            case 3
                % Nakagami-m衰落 (m=2)
                m_nakagami = 2;
                nakagami_temp = sqrt(gamrnd(m_nakagami, 1/m_nakagami, 1, num_samples));
                nakagami_temp = nakagami_temp / sqrt(mean(nakagami_temp.^2)); % 归一化
                temp_results{i} = nakagami_temp;
            case 4
                % Weibull衰落 (k=2)
                k_weibull = 2;
                weibull_temp = wblrnd(1, k_weibull, 1, num_samples);
                weibull_temp = weibull_temp / sqrt(mean(weibull_temp.^2)); % 归一化
                temp_results{i} = weibull_temp;
        end
    end
    
    % 将临时结果复制到最终变量
    rayleigh_samples = temp_results{1};
    rician_samples = temp_results{2};
    nakagami_samples = temp_results{3};
    weibull_samples = temp_results{4};
else
    % 串行计算
    % 瑞利衰落
    rayleigh_samples = sqrt(randn(1,num_samples).^2 + randn(1,num_samples).^2);
    rayleigh_samples = rayleigh_samples / sqrt(mean(rayleigh_samples.^2)); % 归一化

    % 莱斯衰落
    rician_samples = sqrt((sqrt(k_factor_linear) + randn(1,num_samples)).^2 + randn(1,num_samples).^2);
    rician_samples = rician_samples / sqrt(mean(rician_samples.^2)); % 归一化

    % Nakagami-m衰落 (m=2)
    m_nakagami = 2;
    nakagami_samples = sqrt(gamrnd(m_nakagami, 1/m_nakagami, 1, num_samples));
    nakagami_samples = nakagami_samples / sqrt(mean(nakagami_samples.^2)); % 归一化

    % Weibull衰落 (k=2)
    k_weibull = 2;
    weibull_samples = wblrnd(1, k_weibull, 1, num_samples);
    weibull_samples = weibull_samples / sqrt(mean(weibull_samples.^2)); % 归一化
end

%% 统计特性计算
statistical_results = struct();

% 瑞利衰落统计
statistical_results.rayleigh = struct();
statistical_results.rayleigh.mean = mean(rayleigh_samples);
statistical_results.rayleigh.variance = var(rayleigh_samples);
statistical_results.rayleigh.std = std(rayleigh_samples);
statistical_results.rayleigh.envelope = abs(rayleigh_samples);

% 莱斯衰落统计
statistical_results.rician = struct();
statistical_results.rician.mean = mean(rician_samples);
statistical_results.rician.variance = var(rician_samples);
statistical_results.rician.std = std(rician_samples);
statistical_results.rician.envelope = abs(rician_samples);

% Nakagami衰落统计
statistical_results.nakagami = struct();
statistical_results.nakagami.mean = mean(nakagami_samples);
statistical_results.nakagami.variance = var(nakagami_samples);
statistical_results.nakagami.std = std(nakagami_samples);
statistical_results.nakagami.envelope = abs(nakagami_samples);

% Weibull衰落统计
statistical_results.weibull = struct();
statistical_results.weibull.mean = mean(weibull_samples);
statistical_results.weibull.variance = var(weibull_samples);
statistical_results.weibull.std = std(weibull_samples);
statistical_results.weibull.envelope = abs(weibull_samples);

%% 电平通过率 (LCR) 和平均衰落持续时间 (AFD)
thresholds_db = -30:2:0;
thresholds_linear = 10.^(thresholds_db/20);

% 初始化结果
lcr_results = zeros(length(thresholds_db), 4); % 四种衰落类型
afd_results = zeros(length(thresholds_db), 4);

if use_parallel
    % 使用并行计算计算LCR和AFD
    fprintf('使用并行计算计算LCR和AFD...\n');
    
    % 预分配临时结果结构
    temp_results = cell(length(thresholds_linear), 1);
    
    parfor i = 1:length(thresholds_linear)
        level = thresholds_linear(i);
        temp_result = zeros(1, 8); % 4种衰落类型，每种有LCR和AFD两个值
        
        % 瑞利衰落
        rayleigh_envelope = statistical_results.rayleigh.envelope;
        rayleigh_lcr = sum(abs(diff(rayleigh_envelope > level))) / (2*num_samples);
        rayleigh_afd = sum(rayleigh_envelope <= level) / num_samples / (rayleigh_lcr + eps);
        temp_result(1) = rayleigh_lcr;
        temp_result(2) = rayleigh_afd;
        
        % 莱斯衰落
        rician_envelope = statistical_results.rician.envelope;
        rician_lcr = sum(abs(diff(rician_envelope > level))) / (2*num_samples);
        rician_afd = sum(rician_envelope <= level) / num_samples / (rician_lcr + eps);
        temp_result(3) = rician_lcr;
        temp_result(4) = rician_afd;
        
        % Nakagami衰落
        nakagami_envelope = statistical_results.nakagami.envelope;
        nakagami_lcr = sum(abs(diff(nakagami_envelope > level))) / (2*num_samples);
        nakagami_afd = sum(nakagami_envelope <= level) / num_samples / (nakagami_lcr + eps);
        temp_result(5) = nakagami_lcr;
        temp_result(6) = nakagami_afd;
        
        % Weibull衰落
        weibull_envelope = statistical_results.weibull.envelope;
        weibull_lcr = sum(abs(diff(weibull_envelope > level))) / (2*num_samples);
        weibull_afd = sum(weibull_envelope <= level) / num_samples / (weibull_lcr + eps);
        temp_result(7) = weibull_lcr;
        temp_result(8) = weibull_afd;
        
        temp_results{i} = temp_result;
    end
    
    % 将临时结果复制到最终数组
    for i = 1:length(thresholds_linear)
        lcr_results(i, 1) = temp_results{i}(1);  % 瑞利LCR
        afd_results(i, 1) = temp_results{i}(2);  % 瑞利AFD
        lcr_results(i, 2) = temp_results{i}(3);  % 莱斯LCR
        afd_results(i, 2) = temp_results{i}(4);  % 莱斯AFD
        lcr_results(i, 3) = temp_results{i}(5);  % NakagamiLCR
        afd_results(i, 3) = temp_results{i}(6);  % NakagamiAFD
        lcr_results(i, 4) = temp_results{i}(7);  % WeibullLCR
        afd_results(i, 4) = temp_results{i}(8);  % WeibullAFD
    end
else
    % 串行计算
    for i = 1:length(thresholds_linear)
        level = thresholds_linear(i);
        
        % 瑞利衰落
        rayleigh_envelope = statistical_results.rayleigh.envelope;
        rayleigh_lcr = sum(abs(diff(rayleigh_envelope > level))) / (2*num_samples);
        rayleigh_afd = sum(rayleigh_envelope <= level) / num_samples / (rayleigh_lcr + eps);
        lcr_results(i, 1) = rayleigh_lcr;
        afd_results(i, 1) = rayleigh_afd;
        
        % 莱斯衰落
        rician_envelope = statistical_results.rician.envelope;
        rician_lcr = sum(abs(diff(rician_envelope > level))) / (2*num_samples);
        rician_afd = sum(rician_envelope <= level) / num_samples / (rician_lcr + eps);
        lcr_results(i, 2) = rician_lcr;
        afd_results(i, 2) = rician_afd;
        
        % Nakagami衰落
        nakagami_envelope = statistical_results.nakagami.envelope;
        nakagami_lcr = sum(abs(diff(nakagami_envelope > level))) / (2*num_samples);
        nakagami_afd = sum(nakagami_envelope <= level) / num_samples / (nakagami_lcr + eps);
        lcr_results(i, 3) = nakagami_lcr;
        afd_results(i, 3) = nakagami_afd;
        
        % Weibull衰落
        weibull_envelope = statistical_results.weibull.envelope;
        weibull_lcr = sum(abs(diff(weibull_envelope > level))) / (2*num_samples);
        weibull_afd = sum(weibull_envelope <= level) / num_samples / (weibull_lcr + eps);
        lcr_results(i, 4) = weibull_lcr;
        afd_results(i, 4) = weibull_afd;
    end
end

statistical_results.thresholds_db = thresholds_db;
statistical_results.lcr_results = lcr_results;
statistical_results.afd_results = afd_results;

%% 可视化结果
fig_statistics = figure('Name', '统计信道模型分析', 'Position', [250, 250, 1400, 1000]);

% 子图1: 时域波形对比
subplot(3, 3, 1);
time_axis = 0:0.001:0.999;
plot(time_axis, 20*log10(abs(rayleigh_samples(1:1000))), 'LineWidth', 1, 'Color', colors(1, :));
hold on;
plot(time_axis, 20*log10(abs(rician_samples(1:1000))), 'LineWidth', 1, 'Color', colors(2, :));
plot(time_axis, 20*log10(abs(nakagami_samples(1:1000))), 'LineWidth', 1, 'Color', colors(3, :));
plot(time_axis, 20*log10(abs(weibull_samples(1:1000))), 'LineWidth', 1, 'Color', colors(4, :));
grid on;
xlabel('时间 (s)');
ylabel('幅度 (dB)');
title('衰落波形对比');
legend({'瑞利', '莱斯', 'Nakagami', 'Weibull'}, 'Location', 'SouthEast');

% 子图2: 概率密度函数 (PDF)
subplot(3, 3, 2);
[rayleigh_pdf, x_rayleigh] = histcounts(20*log10(abs(rayleigh_samples)), 50, 'Normalization', 'pdf');
[rician_pdf, x_rician] = histcounts(20*log10(abs(rician_samples)), 50, 'Normalization', 'pdf');
[nakagami_pdf, x_nakagami] = histcounts(20*log10(abs(nakagami_samples)), 50, 'Normalization', 'pdf');
[weibull_pdf, x_weibull] = histcounts(20*log10(abs(weibull_samples)), 50, 'Normalization', 'pdf');

plot(x_rayleigh(1:end-1), rayleigh_pdf, 'LineWidth', 2, 'Color', colors(1, :));
hold on;
plot(x_rician(1:end-1), rician_pdf, 'LineWidth', 2, 'Color', colors(2, :));
plot(x_nakagami(1:end-1), nakagami_pdf, 'LineWidth', 2, 'Color', colors(3, :));
plot(x_weibull(1:end-1), weibull_pdf, 'LineWidth', 2, 'Color', colors(4, :));
grid on;
xlabel('幅度 (dB)');
ylabel('概率密度');
title('概率密度函数');
legend({'瑞利', '莱斯', 'Nakagami', 'Weibull'}, 'Location', 'NorthEast');

% 子图3: 累积分布函数 (CDF)
subplot(3, 3, 3);
rayleigh_cdf = cumsum(rayleigh_pdf) * (x_rayleigh(2)-x_rayleigh(1));
rician_cdf = cumsum(rician_pdf) * (x_rician(2)-x_rician(1));
nakagami_cdf = cumsum(nakagami_pdf) * (x_nakagami(2)-x_nakagami(1));
weibull_cdf = cumsum(weibull_pdf) * (x_weibull(2)-x_weibull(1));

plot(x_rayleigh(1:end-1), rayleigh_cdf, 'LineWidth', 2, 'Color', colors(1, :));
hold on;
plot(x_rician(1:end-1), rician_cdf, 'LineWidth', 2, 'Color', colors(2, :));
plot(x_nakagami(1:end-1), nakagami_cdf, 'LineWidth', 2, 'Color', colors(3, :));
plot(x_weibull(1:end-1), weibull_cdf, 'LineWidth', 2, 'Color', colors(4, :));
grid on;
xlabel('幅度 (dB)');
ylabel('累积概率');
title('累积分布函数');
legend({'瑞利', '莱斯', 'Nakagami', 'Weibull'}, 'Location', 'NorthEast');

% 子图4: 电平通过率 (LCR)
subplot(3, 3, 4);
plot(thresholds_db, lcr_results(:, 1), 'LineWidth', 2, 'Color', colors(1, :));
hold on;
plot(thresholds_db, lcr_results(:, 2), 'LineWidth', 2, 'Color', colors(2, :));
plot(thresholds_db, lcr_results(:, 3), 'LineWidth', 2, 'Color', colors(3, :));
plot(thresholds_db, lcr_results(:, 4), 'LineWidth', 2, 'Color', colors(4, :));
grid on;
xlabel('门限电平 (dB)');
ylabel('电平通过率');
title('电平通过率');
legend({'瑞利', '莱斯', 'Nakagami', 'Weibull'}, 'Location', 'NorthEast');

% 子图5: 平均衰落持续时间 (AFD)
subplot(3, 3, 5);
semilogy(thresholds_db, afd_results(:, 1), 'LineWidth', 2, 'Color', colors(1, :));
hold on;
semilogy(thresholds_db, afd_results(:, 2), 'LineWidth', 2, 'Color', colors(2, :));
semilogy(thresholds_db, afd_results(:, 3), 'LineWidth', 2, 'Color', colors(3, :));
semilogy(thresholds_db, afd_results(:, 4), 'LineWidth', 2, 'Color', colors(4, :));
grid on;
xlabel('门限电平 (dB)');
ylabel('平均衰落持续时间');
title('平均衰落持续时间');
legend({'瑞利', '莱斯', 'Nakagami', 'Weibull'}, 'Location', 'NorthEast');

% 子图6: 莱斯因子影响
subplot(3, 3, 6);
k_factors = [-10, -5, 0, 5, 10]; % 不同莱斯因子 (dB)

if use_parallel
    % 并行计算不同莱斯因子的CDF
    fprintf('使用并行计算分析莱斯因子影响...\n');
    parfor i = 1:length(k_factors)
        % 生成不同莱斯因子的样本
        k_linear = 10^(k_factors(i)/10);
        rician_samples_k = sqrt((sqrt(k_linear) + randn(1, num_samples)).^2 + randn(1, num_samples).^2);
        rician_samples_k = rician_samples_k / sqrt(mean(rician_samples_k.^2));
        
        % 计算CDF
        [pdf_k, x_k] = histcounts(20*log10(abs(rician_samples_k)), 50, 'Normalization', 'pdf');
        cdf_k = cumsum(pdf_k) * (x_k(2)-x_k(1));
        
        % 存储结果
        k_factor_results(i).pdf = pdf_k;
        k_factor_results(i).x = x_k;
        k_factor_results(i).cdf = cdf_k;
    end
    
    % 绘制结果
    for i = 1:length(k_factors)
        plot(k_factor_results(i).x(1:end-1), k_factor_results(i).cdf, 'LineWidth', 2, 'Color', colors(i, :));
        hold on;
    end
else
    % 串行计算
    for i = 1:length(k_factors)
        % 生成不同莱斯因子的样本
        k_linear = 10^(k_factors(i)/10);
        rician_samples_k = sqrt((sqrt(k_linear) + randn(1, num_samples)).^2 + randn(1, num_samples).^2);
        rician_samples_k = rician_samples_k / sqrt(mean(rician_samples_k.^2));
        
        % 计算CDF
        [pdf_k, x_k] = histcounts(20*log10(abs(rician_samples_k)), 50, 'Normalization', 'pdf');
    cdf_k = cumsum(pdf_k) * (x_k(2)-x_k(1));
    plot(x_k(1:end-1), cdf_k, 'LineWidth', 2, 'Color', colors(i, :));
        hold on;
    end
end
grid on;
xlabel('幅度 (dB)');
ylabel('累积概率');
title('莱斯因子影响');
legend({sprintf('K=%d dB', k_factors)}, 'Location', 'NorthEast');

% 子图7: 统计特性对比表
subplot(3, 3, 7);
models = {'瑞利', '莱斯', 'Nakagami', 'Weibull'};
means = [statistical_results.rayleigh.mean, statistical_results.rician.mean, ...
         statistical_results.nakagami.mean, statistical_results.weibull.mean];
stds = [statistical_results.rayleigh.std, statistical_results.rician.std, ...
        statistical_results.nakagami.std, statistical_results.weibull.std];

% 创建表格数据
table_data = [means; stds];
imagesc(table_data);
colormap('gray');
colorbar;
set(gca, 'XTick', 1:4, 'XTickLabel', models);
set(gca, 'YTick', 1:2, 'YTickLabel', {'均值', '标准差'});
title('统计特性对比');

% 子图8: 衰落深度对比
subplot(3, 3, 8);
fade_depths = [-3, -6, -10, -15, -20]; % 不同衰落深度 (dB)
fade_probs = zeros(length(fade_depths), 4);

for i = 1:length(fade_depths)
    threshold = 10^(fade_depths(i)/20);
    fade_probs(i, 1) = sum(rayleigh_samples <= threshold) / num_samples;
    fade_probs(i, 2) = sum(rician_samples <= threshold) / num_samples;
    fade_probs(i, 3) = sum(nakagami_samples <= threshold) / num_samples;
    fade_probs(i, 4) = sum(weibull_samples <= threshold) / num_samples;
end

semilogy(fade_depths, fade_probs(:, 1), 's-', 'LineWidth', 2, 'Color', colors(1, :));
hold on;
semilogy(fade_depths, fade_probs(:, 2), 's-', 'LineWidth', 2, 'Color', colors(2, :));
semilogy(fade_depths, fade_probs(:, 3), 's-', 'LineWidth', 2, 'Color', colors(3, :));
semilogy(fade_depths, fade_probs(:, 4), 's-', 'LineWidth', 2, 'Color', colors(4, :));
grid on;
xlabel('衰落深度 (dB)');
ylabel('衰落概率');
title('衰落深度概率');
legend({'瑞利', '莱斯', 'Nakagami', 'Weibull'}, 'Location', 'NorthEast');

% 子图9: 分集增益分析
subplot(3, 3, 9);
% 计算不同分集阶数的性能
snr_range = 0:2:20; % SNR范围 (dB)
ber_rayleigh = 0.5 * (1 - sqrt(snr_range ./ (1 + snr_range))); % 瑞利衰落BER

% 不同分集阶数的BER
orders = [1, 2, 4, 8]; % 分集阶数
for i = 1:length(orders)
    % 简化的分集BER模型
    div_ber = ber_rayleigh.^orders(i);
    semilogy(snr_range, div_ber, 'LineWidth', 2, 'Color', colors(i, :));
    hold on;
end
grid on;
xlabel('SNR (dB)');
ylabel('误码率');
title('分集增益分析');
legend({'L=1', 'L=2', 'L=4', 'L=8'}, 'Location', 'NorthEast');

fprintf('统计信道模型分析完成！\n');

end